#######################################################################################
### Script file to reproduce the analysis in Chapter 4
###
### Created: 6-24-24
### Modified: 
###
#######################################################################################	


## Load the required packages
library(crch)
	library(mnormt)
	library(matrixStats)

## kick out the presidential systems ##
	data <- read.table('Chapter 4/Input/chapter4data.txt', sep = '|', header = TRUE)

## estimate baseline models to create space-holders for bootstrap ##
	
	mod <- crch(sophiaMean ~ lagSophiaMean + cabinet*cmpPoint + coalCab * sophiaPartnersMean + 
		lagCabinet*lagCmpPoint + lagCoalCab * lagSophiaPartners | 
		cabinet + seatShare + cmpSe + ees + cmpChange + log(age),
		data = data)
	summary(mod)

	lin <- lm(sophiaMean ~ lagSophiaMean + cmpPoint + 
		lagCmpPoint, data = data)
	summary(lin)
	
## let's iterate through all of the estimates to recover the full model ##
	set.seed(1)
	data$lagDv <- NA
	data$lagAve <- NA
	data$lagCmpNew <- NA
	data$lagPartnersSophia <- NA
	data$lagPartnersAverage <- NA
	pars <- rmnorm(0, coef(mod), vcov(mod))
	r2 <- c()
	for(i in 1:1000){
		data$dv <- data[, c(12+i)]
		data$cmpNew <- rnorm(data$cmpPoint, data$cmpPoint, data$cmpSe)
		data$partnersSophia <- data[, c(1034+i)]
		
		for(j in 1:length(data$cabinet)){
			if(data$time[j] > 1){
				data$lagDv[j] <- data$dv[data$cmp == data$cmp[j] & data$time == (data$time[j]-1)]
				data$lagAve[j] <- data$average[data$cmp == data$cmp[j] & data$time == (data$time[j]-1)]
				data$lagCmpNew[j] <- data$cmpNew[data$cmp == data$cmp[j] & data$time == (data$time[j]-1)]
				data$lagEes[j] <- data$ees[data$cmp == data$cmp[j] & data$time == (data$time[j]-1)]
				data$lagPartnersSophia[j] <- data$partnersSophia[data$cmp == data$cmp[j] & data$time == (data$time[j]-1)]
				data$lagPartnersAverage[j] <- data$partnersAverage[data$cmp == data$cmp[j] & data$time == (data$time[j]-1)]
				data$lagAge[j] <- data$age[data$cmp == data$cmp[j] & data$time == (data$time[j]-1)]
			}
		}
		
		data$cmpChangeNew <- abs(data$lagCmpNew - data$cmpNew)
		
		mod <- crch(dv ~ lagDv + cabinet*cmpNew + coalCab*partnersSophia + 
			lagCabinet*lagCmpNew + lagCoalCab * lagPartnersSophia | 
			cabinet + seatShare + cmpSe + ees + cmpChangeNew + log(age),
			data = data)

		hold <- rmnorm(100, coef(mod), vcov(mod))
		pars <- rbind(pars, hold)
		r2 <- c(r2, summary(mod)$loglik)
	
		if(match(i, seq(1, 1000, by = 50), nomatch = 0) > 0){
			print(paste('model', i, 'of 1,000 complete...\n'))
		}
	}
	
## write the parameter estimates if you like ##
#	write.table(pars, file = 'chadlPars.txt', sep = '|')
#	pars <- read.table('chadlPars.txt', header = TRUE, sep = '|')
	summary(pars)
	
## output for table 7 ##
	p <- c()
	for(i in 1:dim(pars)[2]){
		if(mean(pars[, i]) > 0){
			p[i] <- length(pars[, i][pars[, i] < 0]) / length(pars[, i])
		}else{
			p[i] <- length(pars[, i][pars[, i] > 0]) / length(pars[, i])
		}
	}

	cbind(round(colMeans(pars), 3), round(colSds(pars), 3), round(p, 3))
	
## output for table 8 ##
## long and short run effects for cabinet parties ##
	steCabCcab <- pars[, 6] + pars[, 12]
	steCabLcab <- pars[, 10] + pars[, 14]
	lteCabCab <- (pars[, 6] + pars[, 10] + pars[, 12] + pars[, 14]) / (1 - pars[, 2])

	steCmpCcab <- pars[, 4] + pars[, 11]
	steCmpLcab <- pars[, 8] + pars[, 13]
	lteCmpCab <- (pars[, 4] + pars[, 8] + pars[, 11] + pars[, 13]) / (1 - pars[, 2])
	mat <- data.frame(steCabCcab, steCabLcab, lteCabCab, steCmpCcab, steCmpLcab, lteCmpCab)
	
	meds <- c()
	lower <- c()
	upper <- c()
	for(i in 1:dim(mat)[2]){
		meds[i] <- median(mat[, i])
		lower[i] <- quantile(mat[, i], 0.025)
		upper[i] <- quantile(mat[, i], 0.975)
	}	

	names <- c('gov STE X', 'gov STE Xt-1', 'gov LTE X', 'cmp STE X', 'cmp STE Xt-1', 'cmp LTE X')
	cbind(names, round(meds, 3), round(lower, 3), round(upper, 3))

## long and short run effects for opposition parties ##
	steCabCopp <- pars[, 6]
	steCabLopp <- pars[, 10]
	lteCabOpp <- (pars[, 6] + pars[, 10]) / (1 - pars[, 2])

	steCmpCopp <- pars[, 4]
	steCmpLopp <- pars[, 8]
	lteCmpOpp <- (pars[, 4] + pars[, 8]) / (1 - pars[, 2])
	mat <- data.frame(steCabCopp, steCabLopp, lteCabOpp, steCmpCopp, steCmpLopp, lteCmpOpp)

	meds <- c()
	lower <- c()
	upper <- c()
	for(i in 1:dim(mat)[2]){
		meds[i] <- median(mat[, i])
		lower[i] <- quantile(mat[, i], 0.025)
		upper[i] <- quantile(mat[, i], 0.975)
	}	

	names <- c('gov STE X', 'gov STE Xt-1', 'gov LTE X', 'cmp STE X', 'cmp STE Xt-1', 'cmp LTE X')
	cbind(names, round(meds, 3), round(lower, 3), round(upper, 3))

## now the variance effects for Figure 8 ##
## first cabinet ##
		
	hold0 <- pars[, 15] + 
		pars[, 16] * 0 + 
		pars[, 17] * mean(data$seatShare, na.rm = TRUE) + 
		pars[, 18] * mean(data$cmpSe, na.rm = TRUE) + 
		pars[, 19] * mean(data$ees, na.rm = TRUE) + 
		pars[, 20] * mean(data$cmpChange, na.rm = TRUE) + 
		pars[, 21] * log(mean(data$age, na.rm = TRUE))
	
	hold1 <- pars[, 15] + 
		pars[, 16] * 1 + 
		pars[, 17] * mean(data$seatShare, na.rm = TRUE) + 
		pars[, 18] * mean(data$cmpSe, na.rm = TRUE) + 
		pars[, 19] * mean(data$ees, na.rm = TRUE) + 
		pars[, 20] * mean(data$cmpChange, na.rm = TRUE) + 
		pars[, 21] * log(mean(data$age, na.rm = TRUE))
	
	cabEffect <- ((exp(hold1) - exp(hold0)) / exp(hold0)) * 100

## ees effect ##
		
	hold0 <- pars[, 15] + 
		pars[, 16] * mean(data$cabinet, na.rm = TRUE) + 
		pars[, 17] * mean(data$seatShare, na.rm = TRUE) + 
		pars[, 18] * mean(data$cmpSe, na.rm = TRUE) + 
		pars[, 19] * 0 + 
		pars[, 20] * mean(data$cmpChange, na.rm = TRUE) + 
		pars[, 21] * log(mean(data$age, na.rm = TRUE))



	hold1 <- pars[, 15] + 
		pars[, 16] * mean(data$cabinet, na.rm = TRUE) + 
		pars[, 17] * mean(data$seatShare, na.rm = TRUE) + 
		pars[, 18] * mean(data$cmpSe, na.rm = TRUE) + 
		pars[, 19] * 1 + 
		pars[, 20] * mean(data$cmpChange, na.rm = TRUE) + 
		pars[, 21] * log(mean(data$age, na.rm = TRUE))
	
	eesEffect <- ((exp(hold1) - exp(hold0)) / exp(hold0)) * 100

## cmp error effect ##
		
	hold0 <- pars[, 15] + 
		pars[, 16] * mean(data$cabinet, na.rm = TRUE) + 
		pars[, 17] * mean(data$seatShare, na.rm = TRUE) + 
		pars[, 18] * mean(data$cmpSe, na.rm = TRUE) + 
		pars[, 19] * mean(data$ees, na.rm = TRUE) +
		pars[, 20] * mean(data$cmpChange, na.rm = TRUE) + 
		pars[, 21] * log(mean(data$age, na.rm = TRUE))
	
	hold1 <- pars[, 15] + 
		pars[, 16] * mean(data$cabinet, na.rm = TRUE) + 
		pars[, 17] * mean(data$seatShare, na.rm = TRUE) + 
		pars[, 18] * (mean(data$cmpSe, na.rm = TRUE) + sd(data$cmpSe, na.rm = TRUE)) + 
		pars[, 19] * mean(data$ees, na.rm = TRUE) +
		pars[, 20] * mean(data$cmpChange, na.rm = TRUE) + 
		pars[, 21] * log(mean(data$age, na.rm = TRUE))
	
	errorEffect <- ((exp(hold1) - exp(hold0)) / exp(hold0)) * 100

## cmp change effect ##
		
	hold0 <- pars[, 15] + 
		pars[, 16] * mean(data$cabinet, na.rm = TRUE) + 
		pars[, 17] * mean(data$seatShare, na.rm = TRUE) + 
		pars[, 18] * mean(data$cmpSe, na.rm = TRUE) + 
		pars[, 19] * mean(data$ees, na.rm = TRUE) +
		pars[, 20] * mean(data$cmpChange, na.rm = TRUE) + 
		pars[, 21] * log(mean(data$age, na.rm = TRUE))
	
	hold1 <- pars[, 15] + 
		pars[, 16] * mean(data$cabinet, na.rm = TRUE) + 
		pars[, 17] * mean(data$seatShare, na.rm = TRUE) + 
		pars[, 18] * mean(data$cmpSe, na.rm = TRUE) + 
		pars[, 19] * mean(data$ees, na.rm = TRUE) +
		pars[, 20] * (mean(data$cmpChange, na.rm = TRUE) + sd(data$cmpChange, na.rm = TRUE)) + 
		pars[, 21] * log(mean(data$age, na.rm = TRUE))
	
	changeEffect <- ((exp(hold1) - exp(hold0)) / exp(hold0)) * 100

## plot the differences ##
	effects <- data.frame(cabEffect, errorEffect, changeEffect, eesEffect)
	
#	pdf('varianceFigure.pdf', width = 8)
	boxplot(effects,
		col = rgb(0, 0, 0, 0.2), 
		border = rgb(0, 0, 0, 1), 
		ylim = c(-30, 40),
		outline = FALSE,
		axes = FALSE,
		xlab = '',
		ylab = 'Percent Change in Perception Variance',
		main = 'Substantive Effects for Perception Variance')
	axis(1, labels = c('Government', 'Manifesto Variance', 'Manifesto Change', 'European Election'),
		at = 1:4)
	axis(2, at = seq(-30, 40, by = 10))
	lines(c(0,5), c(0, 0), lty = 3, col = rgb(0, 0, 0, 0.5))
	text(1, -30, paste('p =', round(1-length(cabEffect[cabEffect < 0])/length(cabEffect), 3)))
	text(2, -30, paste('p =', round(1-length(errorEffect[errorEffect > 0])/length(errorEffect), 3)))
	text(3, -30, paste('p =', round(1-length(changeEffect[changeEffect > 0])/length(changeEffect), 3)))
	text(4, -30, paste('p = 0.000'))
#	dev.off()
	
	